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SUMMARY 

Forced convective diffusion-reaction is considered for viscous axisymmetric extensional 
convecting velocity in the neighborhood of a sphere. For Peclet numbers in the range 0.1 ^ Pe 2 500 
and for Damkohler numbers increasing with increasing Pe but in the overall range 0.02 ^ Da ^ 10, 
average and local Sherwood numbers have been computed. By introducing the eigenfunction expansion 
c(r,0) = £ c n (r)P n (cos0) into the forced convective diffusion equation for the concentration of a 
chemical species undergoing a first order homogeneous reaction and by using properties of the 
Legendre functions P n (cos0), the variable coefficient PDE can be reduced to a system of N+ 1 second 
order ODEs for the radial functions c„ (r), n*0,l,2,...,N. The adaptive grid algorithm of Pereyra and 
Lentini can be used to solve the corresponding 2(N+1) first order differential equations as a two-point 
boundary value problem on 1 i r i r. . Convergence of the expansion for a specific value of N can 
thus be established and provides "spectral" behavior as well as the full concentration field c(r,0). 


INTRODUCTION 

The prevalence of small often spherical or approximately spherical particles, bubbles, or 
droplets in atmospheric physics, chemical reaction engineering, combustion science, and environmental 
technology implies the small Reynolds number (Re « 1) assumed here. For concreteness a solid 
sphere is also assumed. Unlike the axisymmetric uniform streaming motion past a sphere (Stokes, 
1851) that is a reasonable assumption in the neighborhood of sedimenting particles or those in 
fluidized beds, however, the flow field in neighborhood of most particles in other natural, industrial, 
and laboratory circumstances is neither uniform nor can it be assumed to be the so-called slip velocity 
characteristic of the ensemble average over all the particles in complex, even turbulent two phase flow 
such as occurs in stirred tanks, for instance. 

We are interested in considering other physically realistic - and therefore necessarily more 
complicated - flow fields that would have another domain application. The ubiquitous spherical 
geometry and the mathematical simplicity of axisymmetry make the axisymmetric extensional flows 
(Re « 1) a natural candidate. The occurence of extensional flows, in particular of locally 
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axisymmmetric ones in the neighborhood of small spherical particles, bubbles, or drops, one of the 
basic building blocks in the rheology and flow of a wide variety of dispersions. 

There are two axisymmetric extensional flow fields. The biaxial and uniaxial flows both have 
the same streamlines. However, the biaxial flow comes along the axes frotn z - ±« and approaches the 
poles of the sphere symmetrically, departing radially outwardly in the symmetry (x,y)-plane, whereas 
the uniaxial flow is oppositely directed and approaches radially symmetrically in the equatorial plane 
and departs along the ±z axes. Far from the sphere, the dimensionless Cartesian components of the 
velocity are (U „ ,U y ,U , ) - ±(x ,y ,-2z), with ± referring throughout to biaxial and uniaxial, 
respectively. 

For Re-0, all flow fields are at rest, and the Sherwood number is independent of the Peclet 
number and depends solely on the Damkohler number, i.e., Sh-Sh (Da n ). For Re « 1 but not 
identically zero, Sh-Sh (Pe,Da„). Pe no more characterizes convection than Re characterizes the 
velocity field. Different velocity fields convect heat and mass differently, even if they have the same 
small non-zero Re and the same Pe. For Re-0, Sh-l+VDa„, but for Re « 1, although the 
axisymmetric uniform streaming flow and the axisymmetric extensional flows all three have the same 
asymptote for Sh (viz., l+vUa„) as Pe - 0, for Pe « 1 but Pe * 0, the functional dependance upon 
Pe, Da„ will be different for the uniform flow, for the biaxial flow, and for the uniaxial flow, 
Sh=Sh(Pe,Da„) will be different, even though Pe and Da„are identical. What is more, the local mass 
transfer coefficients Sh(0;Pe,Da n ) will be even more different For a uniform streaming flow at 
infinity, Pfeffer and his co-workers have studied homogeneous first order reactions for low Reynolds 
number convective diffusion (Rutland and Pfeffer, 1967), (Chen and Pfeffer, 1970) 

We compare and contrast the results for convective diffusion-reaction for biaxial and uniaxial 
flows with one another and with those for the uniform streaming flow. Our emphasis, however, is on 
the theoretical approach, the mathematical calculations, and the use of the Pereyra-Lentini adaptive 
grid algorithm, above all on certain constraints and computational limitations that arise. 


THEORETICAL APPROACH 

Rather than directly attacking the forced convective diffusioiVdiffusion-reaction equation 
numerically as a variable coefficient partial differential equation in which the extensional velocity field 
introduces the known but complicated set of variable coefficients, we take another tack. We introduce 
the eigenfunction expansion 


c(r, 0) = £ c o^ r ) p n( cos °) 


( 1 ) 


with the P n (cos6) being Legendre functions and the radial functions c„(r) are unknown. By utilizing 
properties of the P n (p), p-cos6, we then reduce the single partial differential equation for c(r,6) to a 
system of N+l ordinary differential equations for the c„(r) and solve them numerically, as outlined in 
the next section. 

The dimensionless forced convective diffusion-reaction equation investigated may be written 

/feU • Vc= V 2 c- Da n c, (2) 

in which the second Damkohler number may be expressed in terms of the first. 


! 
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Da u = Daj Pe , 


(3) 


and the Peclet number Pe for the extensional flow utilizes the characteristic velocity E a , in which E 
is the rate of strain at infinity and a is the radius of the solid sphere: 

Pe - E a 2 / , Da ; - k/ E , Da „ = k a 2 / 2> . (4) 

The low Reynolds number axisymmetric extensional flow has two non-vanishing 
dimensionless velocity components 


U, = ±(r-4r^.|r- )(1 
U fl =±(v-r' 4 )(l - 3 sin 0 cos 

The ± signs refer to the biaxial/uniaxial flows, respectively. The streamlines for the two are identical 
and are shown in Figure 1, with the flow being oppositely directed along the streamlines. The biaxial 
flow comes from infinity toward the poles and exits radially symmetrically in the equatorial plane. The 
axisymmetric extensional creeping flow was obtained by specialization of the solution to the creeping 
flow equation of Batchelor (1970) for a general linear rate of strain at infinity; see also Leal (1992) for 
the final result 

The partial differential equation to be solved, 

Pe (u,(r. 9) + MLlII = v 2 c - Pe Da, c (6) 


- 3 cos 2 9 


(5) 


may be rewritten upon introducing the expansion (1) as 


S' 


dc 


Fir) -3p 2 ) P„{ p) +G(r) c B { r) ( 3p) ( nP„{ p) -nP n _ { ( p) ) 


~P'3F[ r * if) c,(r) 


(7) 


= 0 


in which. 


Fir) .(r - | 

G{r) = (l~r- 5 ). 


r- 2 4 




( 8 ) 
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In order to reduce this to a system of ordinary differential equations by utilizing the 
orthogonality of the P n (|J), we must first reduce all the 0-dependent coefficients to Legendre 
polynomials. To accomplish this we use both algebraic and differential recurrence relations for them 
(Abramowitz and Stegun, 1965) , the former repeatedly as required. Ultimately, the convection terms 
may be written as 


£* 


17*0 


F(r) 


dc_n 

dr 


- 3 F(r) 


dc D 

dr 


- nG(r ) c B ( r ) 


2/7+1 


(/7+l) 2 + n 1 
(2/7+3) (2/7-1) 


- 3 cnrv> a( ' )c ° (r) ) p °^ 
- 3 {*('>-£- -0 ( '> c* r) } ^ M) 

f-lC K) } 


/7(77-l) 

(2/7-1) 


(2/7+1) 


+ C7(r) <?.(/■) 


(9) 


The re main ing terms of the equation need not be rewritten. Upon utilizing the orthogonality of the 
P n (p) and solving for the second derivatives, we obtain for the general n (n * 0,1), 


± Pe 


Fir) dC ° 


dr (2/7+1) 


l) 2 . n 2 ) 

l( 2/7+3) (2/7-1) J 

- 3 (2/1&U h ■ ') Tf* ■ ■< • °- 2) C 1 1 ') '> ) . 

- 3 ( '> - J) Ot ') *.,< r) ) 

g a ( /•) + Zfe,, c a (r) 


( 10 ) 


In the computations and results, it is more informative to vary Pe and Dai (called K in the program 
and figures). 

The boundary conditions on c(r,6) are 


c(r,0) = 1 , 
c(r -*•» , 0 ) = 0 , 


( 11 ) 


which imply 

c 0 (r= 1) = 1, 
c a (r=l) =0, n 2s 1, 
c n { r -•><») =0, n ^ 0. 


( 12 ) 



NUMERICAL ALGORITHM 


The algorithm of Pereyra and Lentini (1978) as codified now in the IMSL subroutine 
DBVPFD was used. It is a robust program for solving two-point boundary value problems. In order to 
solve the ordinary differential equation system represented by (10)-(12), we first must terminate the 
infinite series (1) at N < «•, and the spadal domain at r_ < « . The former leads to a finite system of 
second order equations for which c„(r) s 0 for n < 0, n > N. The latter leads to the modified boundary 
conditions 


c 0 (r=l) = 1 
c n (r = 1) = 0 , n i> 0 
c D (r=r J = 0 , n £ 0 . 


(13) 


I 


I 

I 


The results for c„ (r;Pe,Da,) will obviously depend upon N and r_ . The latter (rj is a parameter that 
can be varied in the program. The former (N) must be selected before the program can be run, but 
once selected (as conservatively as possible), convergence of the series can be established. The other 
crucial computational parameters in the subroutine are the initial and maximum number of mesh points 
(NINIT, MXGRID). 

Finally, the system of N+l second order equations must be converted in the usual way to 
a system of 2(N+1) first order equations in order to employ the IMSL subroutine: 


Co( r ) Ti(*) 


C N - 1 ( r ) ~*y NEgNS , ( x) 
c N (r) ~*y NBjNS ( x) 
dc q dy\ 

= -+ywp& ,U) 

dc i dy r 

1F U) = ^zr (x) 


^N- 1 

dr 


(r) 


dc 


K 


dr 


( r) 


dy neons 

dx 

dy NEONS 


•U) 


dx 


<x) 


~*y NEQNS- 1 ( x ) 
y NEQN^ X) 
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RESULTS AND DISCUSSION 


For a practicing engineer and for many engineering and other scientists and mathematicians, 
the principal goal of such an investigation would be a relation between the average Sherwood number 
(the dimensionless mass transfer coefficient) Sh and the physicochemical parameters, viz., Sh(Pe;Da,). 
Of some practical interest is also the local Sherwood number, which for an axisymmetric converting 
velocity would be expressible as Sh(0;Pe,Da I ), the integral of which, when carried out over the surface 
of the sphere, yields the average Sherwood number Sh. The magnitude of the local Sherwood number 
is the normal derivative of the concentration field c(r,6) at the sphere surface, d c/d r (r,6) | r _, . 
Although the concentration field c(r,0) in other approaches to the forced convective diffusion-reaction 
problem would be the object of the numerical research, it generally receives short schrift as being of 
little practical interest In multipardcle systems, the extent of the concentration fields non-negligible 
level for a single particle can for instance, be useful in assessing, or at least estimating, the minimum 
interparticle distance at which concentration fields of neighboring particles would affect one another. 

We start our discussion, neither with Sh(0;Pe,Daj) nor with c(0;Pe,Da,), but with the object of 
our numerical study, the radial functions c„(r;Pe,Da,), denoted as c„(r;Pe,K). In Figure 2a, b for r. (sR 
in the notation employed throughout the paper) * 10 and Pe*5, K=1 we show the radial functions 
c„(r), for n*0,l,2,...,70 for a biaxial flow. Consistent with the reflection symmetry across the equatorial 
(0*n/2) plane, only the even radial modes are nonvanishing. The radial funrtions decrease in 
magnitude, and N*70 clearly produces a convergent series. 

Biaxial 

When the radial functions are multiplied respectively by their corresponding Legendre 
polynomials, the isocontour plot shown in the upper half of Figure 1 results. The biaxial velocity field 
produces the thin(ner) stagnation concentration boundary layer at the poles. The concentration wake 
then imbeds the equatorial plane symmetrically. There are, to emphasize the point, neither momentum 
boundary layers nor momentum wakes (Re « 1). At the same Pe, r„ , and N, an increase of K from 1 
to 2 reduces (Figure 3) the boundary layer a bit and the wake more, effects that are still more 
pronounced for K*5 (Pe*5) in Figure 4. For K»10 (Pe*5), all of the isocontours (0. 1-0.9 in increments 
of 0.1) except for c*0.01 are spherically symmetric (Figures), as far as is apparent to the naked eye 
(and undoubtedly a boon to theoreticians). 

For an increase of Pe to 50, the K*2 (Figure 6) is of course dissimilar to that for Pe*5, but for 
K*5,10 similar remarks apply to the Pe*50 isocontours: there is one nonspherical isocontour for K=5 
and none at K-10 (Plots not shown). 

For a further increase to Pe*200 (Figure 7) the isocontours show a 2-d salient at K*l, which 
has become almost spherically symmetric at K*2 (Figure 8). For K*5 and 10 (Plots not shown) 
spherical symmetry reigns, the differences being solely the decreasing radii of the circles with 
increasing K. 

The isocontour plot for Pe*500 and K*0.5 (Figure 9) is similar to that for Pe*200 and K-l 
(Figure 7). 

Uniaxial 

The area of stagnation concentration boundary layer for uniaxial extensional creeping flow is 
centered on the stagnation velocity ring, the equator. The concentradon wakes are two, stretching from 
the poles (0*0, n), qualitatively similar to the concentration wake on the downwind pole of a sphere in 
a uniform streaming flow at infinity. Such observations are rendered more faithfully in Figure 10 for 
Pe*5, K*5 than in words. 
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An increase from K*5 to 10 for the same Peclet number (Pe=5) brings about expected 
isocontours (Figure 1 1), as does an increase of Pe to 50 for K=5 (Figure 12) and for K=10 (Figure 
13), by which values spherical isocontours result 

Local Sherwood Numbers 

For fast reactions (Da,*K*5,10), spherical isocontours were observed. An increase in the 
convecdon (i.e., in Pe) served to feed the reaction faster but did not further influence the spherical 
symmetry of the isocontours, once a Pe was reached at which they were spherical. This is nowhere 
more evident than for the biaxial flow in Figures 14a, b for K=5 and 10 respectively; Pe*5, 50, 200, 
500. There are slight local maxima at 0=0, tx and a slight local minimum at 0=n/2. Increases in Pe 
lead to dramatic increases in the level of mass transfer rates without however appreciably affecting 
local values over the surface, relative to one another. The increase in Da, from 5 to 10 increases the 
level of order 10 % for each Pe shown. 

Absent reaction, biaxial convective diffusion produces a local Sherwood number that is peaked 
at 0*0,71 and troughed around 0 *ti/ 2. The clear minimum is reduced rapidly as the maxima increase 
with K (Figure 15a, Pe*5; Figure 15b, Pe*50; Figure 15c, Pe*200; Figure 15d, Pe=500). 

For a uniaxial flow the convective diffusion problem without reaction produces a pronounced 
maximum at 0*tx/2 and minima at 0*0, 7t, as expected (Figure 16a). Also as expected, the strong 
maximum is reduced relative to the minima with increasingly fast reaction (Figure 16a), an effect 
observed with higher Pe (Figure 16b,c). 

Crossplots for K*5 and 10 for the several values of Pe in Figures 17a,b, emphasizing the weak 
0-dependence of Sh(0) for fast reactions. 

Average Sherwood Number 

Different velocity fields convect heat and mass differently, as is evident even for the two types 
of axisymmetric extensional flows. Concentration isocontours, other than those for very high Da, , are 
different for biaxial and uniaxial flows. 

For Pe*5, convective diffusion (K*0 in Tables 1,2) by uniaxial flow manifests a greater 
average mass transfer coefficient than by biaxial flow. Indeed, strictly speaking, for any value of Pe 
and K, Sh(Pe,, Kj) ^ > Sh(Pe,, Kp w , as is evident from Tables 1 and 2. 

Nonetheless, for Pe-5, K*l, Sh^ is greater than Sh* by only 0.07; for K*2, by only 0.03; for 
K*5, by only 0.005. For Pe*50 and K*5, 10, Sh^ > Shy only in the third decimal place, which also 
holds for the same K's , at Pe*200. For K-10, at Pe-500, they differ only in the fourth decimal place. 
Thus, from this limited set of results, Sh is virtually identical for uniaxial and biaxial flows for K*5,10 
for Pea50. For smaller reaction rates and for smaller convection (smaller Pe), small but perceptible 
differences will arise between biaxial and uniaxial creeping flows, with the latter being the larger of 
the two. 
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FIGURE I 


CONCENTRATION ISOCONTOURS 


BIAXIAL FLOW WITH HOMOGENEOUS REACTION 
PE=5 K=1 R=10 ORDER=1 L=70 
LEVELS: 0.01 , 0.1 TO 0.9 BY 0.1 



LEVELS: +/- 0.01, +/-0.1, -5 TO 5 BY 0.5 
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FIGURE 14b 


BIAXIAL FLOW WITH HOMOGENEOUS REACTION 
K-10; PE*5, 50. 200, 500 
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FIGURE 15a 

BIAXIAL FLOW WITH HOMOGENEOUS REACTION 
PE-5; K=0. 5. 10 
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LOCAL SHERWOOD NUMBER 


FIGURE 15b 

BIAXIAL FLOW WITH HOMOGENEOUS REACTION 

PE=50; K= 5, 10 
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FIGURE 15c 

BIAXIAL FLOW WITH HOMOGENEOUS REACTION 

PE=200; K= 5, 10 
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FIGURE 15d 


BIAXIAL FLOW WITH HOMOGENEOUS REACTION 
PE-500; K= 5, 10 



FIGURE 16a 

UNIAXIAL FLOW: HOMOGENEOUS REACTION 
PE«5; K*0, 5, 10 
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LOCAL SHERWOOD NUMBER LOCAL SHERWOOD NUMBER 


FIGURE 16b 

UNIAXIAL FLOW: HOMOGENEOUS REACTION 
PE-50; K= 5, 10 



FIGURE 16c 

UNIAXIAL FLOW: HOMOGENEOUS REACTION 
PE-200; K- 5, 10 
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FIGURE 17a 


UNIAXIAL FLOW: HOMOGENEOUS REACTION 
K>5: PE»5, 50, 200 
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FIGURE 17b 


UNIAXIAL FLOW: HOMOGENEOUS REACTION 
K=10; PE»5, 50. 200, 500 
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TABLE 1: Average Sherwood Numbers for Biaxial Flow 


Pe 

K 

R 

Avg. Sherwood N 

NI 

NM 

NF 

NE 

N 

TOL 

5 

0 

10 

2.4222332546 

60 

600 

166 

ca 

70 

IE-06 

5 

1 

10 

3.4844020321 

60 


135 

Rl 

70 

IE-06 

5 

2 

10 

4.2826214867 

60 


135 

142 

70 

IE-06 

5 

5 

10 

6.0326632737 

60 


104 

142 

70 

IE-06 

5 

10 

10 

8.0806645862 

60 


101 

142 

70 

IE-06 

50 

5 

5 

16.8417150406 

100 

875 

140 

142 

70 

IE-06 

50 

10 

5 

23.3676877824 

100 

875 

130 

142 

70 

IE-06 

200 

5 

5 

32.6445658252 

100 

875 

219 

142 

70 

IE-06 

200 

10 

5 

45.7258627584 

100 

875 

235 

142 

70 

IE-06 

500 

5 

5 

51.0160617125 

100 

875 

425 

142 

70 

IE-06 

500 

10 

5 

71.7138264523 

100 

875 

202 

142 

70 

IE-06 


TABLE 2: Average Sherwood Numbers for Uniaxial Flow 


Pe 

K 

R 

Avg. Sherwood N 

NI 

NM 

NF 

NE 

N 

TOL 

5 

0 

10 

2.6345022231 

60 

600 

141 

142 

70 

IE-06 

5 

1 

10 

3.5533852471 

60 


178 

142 

70 

IE-06 

5 

2 

10 

4.3116939874 

60 


189 

142 

70 

IE-06 

5 

5 

10 

6.0374513084 

60 


154 

142 

70 

IE-06 

5 

10 

10 

8.0814290608 

60 


141 

142 

70 

IE-06 

50 

5 


16.8450206823 

100 

875 

309 

142 

70 

IE-06 

50 

10 


23.3680425224 

100 

875 

243 

142 

70 

IE-06 

200 

5 


32.6462378512 

100 

875 

527 

142 

70 

IE-06 

200 

10 


45.7260072732 

100 

875 

353 

142 

70 

IE-06 

500 

10 


71.7138969816 

100 

875 

417 

142 

70 

IE-06 


NI Number of initial grid points, including the endpoints (NINIT) 
NM Maximum number of grid points allowed (MXGRID) 

NF Number of final grid points, including the endpoints (NFINAL) 
NE Number of (first order) differential equations (NEQNS) 

N Number of terms in the cigenfimction expansion 

TOL Relative error control parameter 
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